#This program processes soil temperature date given on http://www.ncdc.noaa.gov/oa/soil/data/ in order
#to obtain the following: average annual soil temperature at depth 50 cm from the soil surface.
#To obtain this average, daily values are summed to give a monthly average. Then monthly averages
#are summed to give annual averages.
#Final output files will have columns of station id, latitude, longitude, year, and annual soil temperature value. 
#
#Program by H.Edwin Winzeler
#Completed December 12, 2011
import sys, os
import math
from datetime import datetime, timedelta
from string import *
NoDataVal = 999
#check for proper number of command line arguments
if (len(sys.argv)!= 4):
    sys.stderr.write ("Usage: %s <StateIdfile> <StationDataFile> <Outfile>.\n\n" % sys.argv[0])
    sys.exit()
#handle command line arguments
StnFile = sys.argv[1]
DataFile = sys.argv[2]
OutFile = sys.argv[3]
if os.path.isfile ( StnFile ):#obtain the station id data and its contents
    #Station ID file is found, so open it and read its contents
    fin = open( StnFile, "r" )
    lines = fin.readlines()
    fin.close()
StLN = []#Initialize the station line length list. It shows how many lines of data each station has.
#parse the station id file and populate the station line length list.
for lidx in range(6,len(lines)):
    lines[lidx] = lines[lidx].strip().split()
    lines[lidx][0] = (lines[lidx][0])
    lines[lidx][3] = int(lines[lidx][3])
    lines[lidx][4] = float(lines[lidx][4])
    lines[lidx][5] = (lines[lidx][5]).strip('N')#Remove unnecessary northing symbol
    lines[lidx][5] = float(lines[lidx][5])#Convert value to float
    lines[lidx][5] = (lines[lidx][5])/60#Convert the latitude minutes to decimal degrees
    lines[lidx][4] = lines[lidx][4]+lines[lidx][5]#Add the latitude minutes to the lat degrees
    lines[lidx][6] = float(lines[lidx][6])
    lines[lidx][7] = (lines[lidx][7]).strip('W')#Remove unnecessary westing symbol
    lines[lidx][7] = float(lines[lidx][7])#Convert to float
    lines[lidx][7] = (lines[lidx][7])/60#Convert longitude minutes to decimal degrees
    lines[lidx][6] = lines[lidx][6]+lines[lidx][7]#Add the longitude minutes to the longitude degrees
    StLN.append(lines[lidx][3])#Add station line lengths to the station line-length list
StationInfo = [[0]]*len(StLN)
#POPULATE STATION LINE DATA. THIS GIVES INFORMATION ABOUT EACH STATION, ITS GEOGRAPHIC LOCATION AND NUMBER OF RECORDS.
for lidx in range(len(StLN)):
    b = (lidx + 6)
    StationInfo[lidx] = [999]*3
    StationInfo[lidx][0] = lines[b][0]
    StationInfo[lidx][1] = lines[b][4]
    StationInfo[lidx][2] = lines[b][6]
#OPEN THE DATA FILE AND READ IT
if os.path.isfile ( DataFile ):
    #Station ID file is found, so open it and read its contents
    fin = open( DataFile, "r" )
    dlines = map( rstrip, fin.readlines())
    fin.close()
#DELETE REDUNDANT STATION DATA. SINCE THE STATION DATA LINES HAVE THE FIRST LINE EQUAL TO THE STATION DATA ID, WE DON'T WANT TO PARSE THAT LINE, BUT THE OTHER LINES ACCORDING TO THEIR POSITION IN THE DLINES AS GIVEN IN StLN matrix. I will delete the lines that give station data ID, since I already have it in the 
DelLine = 0
Initial_dlines = len(dlines)
for dlidx in range(len(StLN)): #This loop deletes all the lines that contain
    #station ID information. Cleans up the station data so that it only
    #contains data, and no station information. We have the station id information
    #in the StationInfo array, so we don't need it to be duplicated.
    del dlines[DelLine]        
    DelLine = DelLine+(StLN[dlidx])
Final_dlines = len(dlines)
NumDeletedLines = (Initial_dlines)-(Final_dlines)
Data = [0]*len(dlines)#INITIALIZE THE DATA ARRAY
Flags = [0]*len(dlines)
#POPULATE THE DATA ARRAY WITH DATA FROM THE INPUT DATA FILE
for lidx in range(len(dlines)):
    Data[lidx]=['DataArb']*33# 'DataArb' is an arbitrary string useful for identifying gaps in the data array.
    Flags[lidx]=['flags']*33 #The value 'flags' is an arbitrary string so I can identify areas where the array has not filled in with data.
    Year = int(dlines[lidx][0:4])
    Month = int(dlines[lidx][4:6])
    Data[lidx][0]=datetime(Year,Month,1)
    Data[lidx][1]=int(dlines[lidx][7:8])#This parses the soil depth code. We need data lines that have the
    #soil depth code 4, which indicates measurements made at 50 cm. All other observations will not be used
    #in the annual averages.
    Flags[lidx][0]=datetime(Year,Month,1)
    Flags[lidx][1]=int(dlines[lidx][7:8])
    startpos = 11
    for idx in range(31):#31 days in each month, padded with  non-data values (999)
        endpos = startpos + 6 #the end position of the data phrase, which consists of a temperature reading and a flag value
        Temposend = startpos + 3#the end position of the temperature part of the phrase
        flagposbeg = Temposend #the beginning position of the data flag
        flagposend = Temposend+3 #the end position of the data flag
        Data[lidx][2+idx] = int(dlines[lidx][startpos:Temposend])#Even though I am going to get the average value per month and per year, there is no need for floating point precision because we only have two sig digits to start with, so adding precision is illusory. (No need for float)
        Flags[lidx][2+idx] = int(dlines[lidx][flagposbeg:flagposend])
        startpos=endpos
#GET RID OF PROBLEMATIC DATA, REPLACE IT WITH NON-DATA 999 VALUES
for lidx in range(len(dlines)): #This changes all data that has problem flags to non-data (that is, data that has internal inconsistencies,
    #failures of the spice check, failures of the limited check, failures of the rate-of-chance limeit, or missing values). See dataparameters-2.doc on http://www.ncdc.noaa.gov/oa/soil/data/ for details.
    NumData = len(Data[lidx])
    for ix in range(NumData-2):#Because we don't want to overite the 1 and 2nd compenents of the data line (these are the date and the depth codes)
        if Flags[lidx][ix+2] != 0:#Because we don't want to overite the 1 and 2nd compenents of the data line (these are the date and the depth codes)
            Data[lidx][ix+2] = 999
#Get the mean value for the Temperature for the month and ignore non-data values
def get_mean(values, N = "", NoData=NoDataVal, Skip = ""): #This function was taken from the program stats.py by Keith Cherkauer
	if not N: N = len(values)
	mean = 0
	Nact = 0
	Nskip = 0
	for i in range(N):
		if values[i] != NoData:
		    if Skip and values[i] == Skip:
			Nskip = Nskip +1
		    else:
		        mean = mean + values[i]
		    Nact = Nact +1
	if Nact-Nskip>0:
	    mean = mean / (Nact - Nskip)
	else:
	    mean = NoData
	return mean, Nact
#GET THE MEAN MONTHLY TEMPERATURE FOR EACH DATA LINE
DataNoDates = [999]*len(Data)
MeanMonthly = [999]*len(Data)
ValidMeasurements = [999]*len(Data)
for lidx in range(len(Data)):
    DataNoDates=Data[lidx][2:33]
    get_mean(DataNoDates)
    MeanMonthly[lidx]=[999]*3
    MeanMonthly[lidx][0]=Data[lidx][0]
    MeanMonthly[lidx][1]=Data[lidx][1]
    MeanMonthly[lidx][2]=get_mean(DataNoDates)[0]
    ValidMeasurements[lidx]=[999]*3
    ValidMeasurements[lidx][0] = Data[lidx][0]
    ValidMeasurements[lidx][1] = Data[lidx][1]
    ValidMeasurements[lidx][2] = get_mean(DataNoDates)[1]

#This gives us measurements at all depths, but we only need those for the 50 cm depth
MeanAnnual = [999]*len(Data)
NumYears = (Data[-1][0].year) - (Data[0][0].year)

#I only want soil moisture observations that occur at 50 cm, or depth 4 as listed in column 4 (which is Data[lidx][1]) 
#COUNT HOW MANY STATIONS THERE ARE IN THE GIVEN STATE THAT HAVE SOIL TEMPERATURE DATA AT 50 CM
ObsAt50 = 0
for lidx in range(len(Data)):
    if Data[lidx][1] == 4:
        ObsAt50+=1
Obs50Lines = []#INITIALIZE THE ARRAY, ASSUMING THAT ALL LINES GIVE DATA AT 50 CM. LATER WE WILL DELETE THE EXTRA LINES.
#POPULATE THE Obs50Lines WITH ONLY DATA FROM THE 50 CM DEPTH.
for lidx in range(len(Data)):
    if Data[lidx][1] == 4:
        Obs50Lines.append(Data[lidx])
        ObsAt50+=1
print "Obs50Lines,", Obs50Lines
print "The number of monthly observations at depth 50 cm is,", ObsAt50

#GET THE MEAN MONTHLY TEMPERATURE FOR EACH DATA LINE AT 50 CM
DataNoDates50 = [999]*len(Obs50Lines)
MeanMonthly50 = [999]*len(Obs50Lines)
ValidMeasurements50 = [999]*len(Obs50Lines)
for lidx in range(len(Obs50Lines)):
    DataNoDates50=Obs50Lines[lidx][2:33]
    get_mean(DataNoDates50)
    MeanMonthly50[lidx]=[999]*3
    MeanMonthly50[lidx][0]=Obs50Lines[lidx][0]
    MeanMonthly50[lidx][1]=Obs50Lines[lidx][1]
    MeanMonthly50[lidx][2]=get_mean(DataNoDates50)[0]
    ValidMeasurements50[lidx]=[999]*3
    ValidMeasurements50[lidx][0] = Obs50Lines[lidx][0]
    ValidMeasurements50[lidx][1] = Obs50Lines[lidx][1]
    ValidMeasurements50[lidx][2] = get_mean(DataNoDates50)[1]

print "MeanMonthly50,",MeanMonthly50

#WHERE THERE ARE MORE THAN ONE OBSERVATION FOR A GIVEN MONTH AT A GIVEN STATION, AVERAGE THEM. LATER WE AVERAGE ALL MONTHS TOGETHER FOR EACH YEAR
#TO GIVE US AN ANNUAL SOIL TEMPERATURE VALUE.

##print "NumDeletedLines,", NumDeletedLines
##print "len(dlines),", len(dlines)
##print "dlines[0:1],", dlines[0:1]
StObs = [999]*(len(lines)-6)#initialize the station array
for lidx in range(len(lines)-6):
    Columns = 4
    StObs[lidx]=[999]*Columns
    StObs[lidx][0]=lines[lidx+6][0]
    StObs[lidx][1]=lines[lidx+6][3]
    StObs[lidx][2]=lines[lidx+6][4]
    StObs[lidx][3]=lines[lidx+6][6]
##print "StObs",StObs
    
fout = open(OutFile, "w")
fout.write("Station\tMonthsOfObservation\tLatitude\tLongitude\t\n")
for lidx in range(len(StObs)):
    fout.write("%s" % StObs[lidx])
    fout.write("\n")
